home *** CD-ROM | disk | FTP | other *** search
/ Collection of Tools & Utilities / Collection of Tools and Utilities.iso / ada / gnat1792.zip / gnat179b / t-adainc / a-ngcefu.adb < prev    next >
Text File  |  1994-05-19  |  18KB  |  669 lines
  1. ------------------------------------------------------------------------------
  2. --                                                                          --
  3. --                         GNAT RUNTIME COMPONENTS                          --
  4. --                                                                          --
  5. --                     A D A . N U M E R I C S . G C E F                    --
  6. --                                                                          --
  7. --                                 B o d y                                  --
  8. --                                                                          --
  9. --                            $Revision: 1.1 $                              --
  10. --                                                                          --
  11. --           Copyright (c) 1992,1993,1994 NYU, All Rights Reserved          --
  12. --                                                                          --
  13. -- GNAT is free software;  you can  redistribute it  and/or modify it under --
  14. -- terms of the  GNU General Public License as published  by the Free Soft- --
  15. -- ware  Foundation;  either version 2,  or (at your option) any later ver- --
  16. -- sion.  GNAT is distributed in the hope that it will be useful, but WITH- --
  17. -- OUT ANY WARRANTY;  without even the  implied warranty of MERCHANTABILITY --
  18. -- or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License --
  19. -- for  more details.  You should have  received  a copy of the GNU General --
  20. -- Public License  distributed with GNAT;  see file COPYING.  If not, write --
  21. -- to the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. --
  22. --                                                                          --
  23. ------------------------------------------------------------------------------
  24.  
  25. --  This is an early trial implementation, simplified for early gnat.
  26. --  It is not a "strict" implementation.
  27. --  All Ada required exception handling is provided.
  28. --  Many special cases are handled locally to avoid unnecessary calls
  29.  
  30. with Ada.Numerics.Generic_Elementary_Functions;
  31.  
  32. package body Ada.Numerics.Generic_Complex_Elementary_Functions is
  33.  
  34.    package Elementary_Functions is new
  35.       Ada.Numerics.Generic_Elementary_Functions (Real'Base);
  36.    use Elementary_Functions;
  37.  
  38.    PI : constant := 3.14159_26535_89793_23846_26433_83279_50288_41971;
  39.    PI_2 : constant := PI / 2.0;
  40.    Log_Two : constant := 0.69314_71805_59945_30941_72321_21458_17656_80755;
  41.  
  42.    Epsilon : Real'Base;
  43.    Square_Root_Epsilon     : Real'Base;
  44.    Inv_Square_Root_Epsilon : Real'Base;
  45.    Root_Root_Epsilon       : Real'Base;
  46.    Log_Inverse_Epsilon_2   : Real'Base;
  47.  
  48.    Complex_Zero : constant Complex := Compose_From_Cartesian (0.0,  0.0);
  49.    Complex_One  : constant Complex := Compose_From_Cartesian (1.0,  0.0);
  50.    Complex_I    : constant Complex := Compose_From_Cartesian (0.0,  1.0);
  51.    HALF_PI      : constant Complex := Compose_From_Cartesian (PI_2, 0.0);
  52.  
  53.    ----------
  54.    -- Sqrt --
  55.    ----------
  56.  
  57.    function Sqrt (X : Complex) return Complex is
  58.       Z   : Complex := X; -- remove if definition gets fixed
  59.       R   : Real'Base;
  60.       R_X : Real'Base;
  61.       R_Y : Real'Base;
  62.       XR  : Real'Base := abs Re (Z);
  63.       YR  : Real'Base := abs Im (Z);
  64.       A   : Real'Base;
  65.    begin
  66.  
  67.       if XR > YR then
  68.          A := YR / XR;
  69.          A := A * A;
  70.  
  71.          if A < 32.0 * Square_Root_Epsilon then
  72.  
  73.             if Re (Z) > 0.0 then
  74.                R_X := Sqrt (XR + 0.5  *
  75.                            (0.5 * YR * YR / XR - (YR * YR / XR) * A / 8.0));
  76.                R_Y := Sqrt (0.5 *
  77.                            (0.5 * YR * YR / XR - (YR * YR / XR) * A / 8.0));
  78.             else
  79.                R_X := Sqrt (0.5 *
  80.                            (0.5 * YR * YR / XR - (YR * YR / XR) * A / 8.0));
  81.                R_Y := Sqrt (XR + 0.5 *
  82.                            (0.5 * YR * YR / XR - (YR * YR / XR) * A / 8.0));
  83.             end if;
  84.  
  85.          else
  86.             R  := XR * Sqrt (1.0 + A);
  87.             R_X := Sqrt (0.5 *  (R + Re (Z)));
  88.             R_Y := Sqrt (0.5 *  (R - Re (Z)));
  89.          end if;
  90.  
  91.       else                                      -- YR > XR
  92.          A := XR / YR;
  93.          A := A * A;
  94.  
  95.          if A < 32.0 * Square_Root_Epsilon then -- YR >> XR
  96.             R  := YR + 0.5 * XR * XR / YR - 0.125 *  (XR * XR / YR) * A;
  97.             R_X := Sqrt (0.5 *  (YR + (0.5 * XR * XR / YR + Re (Z))));
  98.             R_Y := Sqrt (0.5 *  (YR + (0.5 * XR * XR / YR - Re (Z))));
  99.          else
  100.             R  := YR * Sqrt (1.0 + A);
  101.             R_X := Sqrt (0.5 *  (R + Re (Z)));
  102.             R_Y := Sqrt (0.5 *  (R - Re (Z)));
  103.          end if;
  104.       end if;
  105.  
  106.       if Im (Z) < 0.0 then                 -- halve angle, Sqrt of magnitude
  107.          R_Y := -R_Y;
  108.       end if;
  109.       return Compose_From_Cartesian (R_X, R_Y);
  110.  
  111.    exception
  112.       when Constraint_Error =>
  113.          R := Modulus (Compose_From_Cartesian (Re (Z / 4.0), Im (Z / 4.0)));
  114.          R_X := 2.0 * Sqrt (0.5 * R + 0.5 * Re (Z / 4.0));
  115.          R_Y := 2.0 * Sqrt (0.5 * R - 0.5 * Re (Z / 4.0));
  116.  
  117.          if Im (Z) < 0.0 then -- halve angle, Sqrt of magnitude
  118.             R_Y := -R_Y;
  119.          end if;
  120.  
  121.          return Compose_From_Cartesian (R_X, R_Y);
  122.    end Sqrt;
  123.  
  124.    ---------
  125.    -- Log --
  126.    ---------
  127.  
  128.    function Log (X : Complex) return Complex is
  129.       Z : Complex := X;
  130.       RE_Z, IM_Z : Real'Base;
  131.    begin
  132.  
  133.       if Re (Z) = 0.0 and then Im (Z) = 0.0 then
  134.          raise Constraint_Error;
  135.  
  136.       elsif abs (1.0 - Re (Z)) < Root_Root_Epsilon and then
  137.                       abs Im (Z) < Root_Root_Epsilon then
  138.          Set_Re (Z, Re (Z) - 1.0);
  139.          return (1.0 - (1.0 / 2.0 -
  140.                        (1.0 / 3.0 - (1.0 / 4.0) * Z) * Z) * Z) * Z;
  141.       end if;
  142.  
  143.       begin
  144.          RE_Z := Log (Modulus (Z));
  145.       exception
  146.          when Constraint_Error =>
  147.             RE_Z := Log (Modulus (Z / 2.0)) - Log_Two;
  148.       end;
  149.  
  150.       IM_Z := Arctan (Im (Z), Re (Z));
  151.  
  152.       if IM_Z > PI then
  153.          IM_Z := IM_Z - 2.0 * PI;
  154.       end if;
  155.  
  156.       return Compose_From_Cartesian (RE_Z, IM_Z);
  157.    end Log;
  158.  
  159.    ---------
  160.    -- Exp --
  161.    ---------
  162.  
  163.    function Exp (X : Complex) return Complex is
  164.       Z : Complex := X;
  165.       EXP_RE_Z : Real'Base := Exp (Re (Z));
  166.    begin
  167.       return Compose_From_Cartesian (EXP_RE_Z * Cos (Im (Z)),
  168.                                      EXP_RE_Z * Sin (Im (Z)));
  169.    end Exp;
  170.  
  171.  
  172.    function Exp (X : Imaginary) return Complex is
  173.       IM_Z : Real'Base := Im (X);
  174.    begin
  175.       return Compose_From_Cartesian (Cos (IM_Z), Sin (IM_Z));
  176.    end Exp;
  177.  
  178.    --------
  179.    -- ** --
  180.    --------
  181.   
  182.    function "**"
  183.      (Left : Complex;
  184.      Right : Complex)
  185.    return Complex is
  186.       Z1 : Complex := Left;
  187.       Z2 : Complex := Right;
  188.    begin
  189.       if Re (Z2) = 0.0 and then Im (Z2) = 0.0 and then
  190.          Re (Z1) = 0.0 and then Im (Z1) = 0.0 then
  191.          raise Argument_Error;
  192.  
  193.       elsif Re (Z1) = 0.0 and then Im (Z1) = 0.0 and then
  194.              Re (Z2) < 0.0 then
  195.          raise Constraint_Error;
  196.  
  197.       elsif Re (Z1) = 0.0 and then Im (Z1) = 0.0 then
  198.          return Z1;
  199.  
  200.       elsif Re (Z2) = 0.0 and then Im (Z2) = 0.0 then
  201.          return 1.0 + Z2;
  202.  
  203.       elsif Re (Z2) = 1.0 and then Im (Z2) = 0.0 then
  204.          return Z1;
  205.  
  206.       else
  207.          return Exp (Z2 * Log (Z1));
  208.       end if;
  209.    end "**";
  210.  
  211.  
  212.    function "**"
  213.      (Left : Real'Base;
  214.       Right : Complex)
  215.    return Complex is
  216.       X : Real'Base := Left;
  217.       Z : Complex := Right;
  218.    begin
  219.       if Re (Z) = 0.0 and then Im (Z) = 0.0 and then
  220.          X = 0.0 then
  221.          raise Argument_Error;
  222.  
  223.       elsif X = 0.0 and then Re (Z) < 0.0 then
  224.          raise Constraint_Error;
  225.  
  226.       elsif X = 0.0 then
  227.          return Compose_From_Cartesian (X, 0.0);
  228.  
  229.       elsif Re (Z) = 0.0 and then Im (Z) = 0.0 then
  230.          return Complex_One;
  231.  
  232.       elsif Re (Z) = 1.0 and then Im (Z) = 0.0 then
  233.          return Compose_From_Cartesian (X, 0.0);
  234.  
  235.       else
  236.          return Exp (Log (X) * Z);
  237.       end if;
  238.    end "**";
  239.  
  240.  
  241.    function "**" (Left : Complex;
  242.                    Right : Real'Base) return Complex is
  243.       Z : Complex := Left;
  244.       Y : Real'Base := Right;
  245.    begin
  246.       if  Y = 0.0 and then
  247.           Re (Z) =